## Perspectives on Politics
## Rally for Democracy Script 

## Have reordered this file as figures and tables appear in the final manuscript. 


library(ggplot2)  
library(scales)  
library(ggrepel)
library(WDI)
library(data.table)
library(plyr)
library(readstata13)  
library(gridExtra)
library(gnumeric)
library(DataCombine) 
library(arm)
library(haven)
library(stargazer)
library(lme4)
library(lmerTest)



theme_article <- function(base_size=11, font="Times New Roman"){
  
  txt <- element_text(size = base_size+2, colour = "black", face = "plain", family=font)
  bold_txt <- element_text(size = base_size+2, colour = "black", face = "bold", family=font)
  
  theme_bw(base_size = base_size, base_family = font)+
    theme(
      legend.key = element_blank(), 
      strip.background = element_blank(), 
      text = txt, 
      plot.title = txt, 
      axis.title = txt, 
      axis.text = txt, 
      legend.title = bold_txt, 
      legend.text = txt ) +
    
    ############## lines, background, panel
    theme(
      panel.grid.major.x = element_blank() , 
      panel.grid.minor.x = element_blank(), 
      panel.grid.major.y = element_line( size=.5, color="#CCCCCC"), panel.grid.minor.y = element_line( size=.5, color="#DDDDDD"),
      
      axis.line.y = element_line(colour = "#CCCCCC", size = 0.3),
      axis.line.x = element_line(colour = "#CCCCCC", size = 0.3),
      #### remove Tick marks
      axis.ticks=element_blank(),
      
      ### no legend!
      legend.position = "none", 
      
      ## background
      plot.background = element_rect(fill = "white",colour = "grey",size = 0.5)
    )
}


# a fork of ggsave: not that important but is the function invoked in this script. 
save.plot <- function(name = "test.png", width=6, height= 6, usewd=F){
  wd <- getwd() 
  if(usewd == F) {
    setwd("~/R") }
  else {setwd(wd)}
  ggsave(name, width= width, height= height) 
  dev.off()
  setwd(wd) # reset the working directory to what it was before 
}

# writing a quick save.tiff function for scalable print versions 
save.tiff <- function(plot, name = "test.tiff", resnum=300, w=6, h=6, usewd=F){
  wd <- getwd()
  if(usewd == F) {
    setwd("~/R") }
  else {setwd(wd)}
tiff(filename = name, res = resnum, width = w, height = h, units = 'in', pointsize = 2)
plot(plot)
dev.off() 
setwd(wd)
}


## shorthand for as.data.frame
adf <- function(x) { 
  return(as.data.frame(x)) 
}


# bespoke sort function 
stort <- function(dataset,sortover) {
  return ( dataset[order(sortover) , ] )
} 


## DATA 

# this is the new file that includes 2023 Pew data as well. 
# we use this version of the file for most updated final models and new Appendix figures.

pewdata <- read_dta("Pew2010-2023_VDEM.dta")


# make an age square variable 
pewdata$age_square <- pewdata$age * pewdata$age 

# these still needed recoding for polarity 
pewdata$fav_nato <- 5 - pewdata$fav_nato
pewdata$fav_un <- 5 - pewdata$fav_un



## ideology variable.

pewdata$extremist <- NA
pewdata$extremist[pewdata$ideology==0] <- 1
pewdata$extremist[pewdata$ideology==6] <- 1
pewdata$extremist[pewdata$ideology==1] <- 0
pewdata$extremist[pewdata$ideology==2] <- 0
pewdata$extremist[pewdata$ideology==3] <- 0
pewdata$extremist[pewdata$ideology==4] <- 0
pewdata$extremist[pewdata$ideology==5] <- 0

pewdata$leftist <- pewdata$extremist 
pewdata$leftist[pewdata$ideology==0] <- 0

pewdata$rightist <- pewdata$extremist 
pewdata$rightist[pewdata$ideology==6] <- 0


pewdata$centrist <- NA
pewdata$centrist[pewdata$ideology==0] <- 0
pewdata$centrist[pewdata$ideology==6] <- 0
pewdata$centrist[pewdata$ideology==1] <- 0
pewdata$centrist[pewdata$ideology==2] <- 1
pewdata$centrist[pewdata$ideology==3] <- 1
pewdata$centrist[pewdata$ideology==4] <- 1
pewdata$centrist[pewdata$ideology==5] <- 0



pewdata$globalsouth <- NA 

# few countries were commented out to leave as NA due to low sample inclusion.
# they would be excluded from sample in final figures, 
# though, as eventually we added a constant country variable to be sure: 
# this step became redundant, but you can compare with both. 

pewdata$globalsouth[pewdata$country=="Argentina"] <- 1
pewdata$globalsouth[pewdata$country=="Australia"] <- 0
# pewdata$globalsouth[pewdata$country=="Bangladesh"] <- 1
pewdata$globalsouth[pewdata$country=="Belgium"] <- 0
pewdata$globalsouth[pewdata$country=="Bolivia"] <- 1
pewdata$globalsouth[pewdata$country=="Brazil"] <- 1
pewdata$globalsouth[pewdata$country=="Britain"] <- 0
pewdata$globalsouth[pewdata$country=="Bulgaria"] <- 0
# pewdata$globalsouth[pewdata$country=="Burkina Faso"] <- 1
pewdata$globalsouth[pewdata$country=="Canada"] <- 0
pewdata$globalsouth[pewdata$country=="Chile"] <- 1
# pewdata$globalsouth[pewdata$country=="China"] <- 1
# pewdata$globalsouth[pewdata$country=="Colombia"] <- 1
pewdata$globalsouth[pewdata$country=="Czech Republic"] <- 0
# pewdata$globalsouth[pewdata$country=="Denmark"] <- 0
pewdata$globalsouth[pewdata$country=="Egypt"] <- 1
pewdata$globalsouth[pewdata$country=="El Salvador"] <- 1
# pewdata$globalsouth[pewdata$country=="Ethiopia"] <- 1
pewdata$globalsouth[pewdata$country=="France"] <- 0
pewdata$globalsouth[pewdata$country=="Germany"] <- 0
pewdata$globalsouth[pewdata$country=="Ghana"] <- 1
pewdata$globalsouth[pewdata$country=="Greece"] <- 0
pewdata$globalsouth[pewdata$country=="Hungary"] <- 0
pewdata$globalsouth[pewdata$country=="India"] <- 1
pewdata$globalsouth[pewdata$country=="Indonesia"] <- 1
pewdata$globalsouth[pewdata$country=="Israel"] <- 0
pewdata$globalsouth[pewdata$country=="Italy"] <- 0
pewdata$globalsouth[pewdata$country=="Japan"] <- 0
pewdata$globalsouth[pewdata$country=="Jordan"] <- 1
pewdata$globalsouth[pewdata$country=="Kenya"] <- 1
pewdata$globalsouth[pewdata$country=="Lebanon"] <- 1
pewdata$globalsouth[pewdata$country=="Lithuania"] <- 0
pewdata$globalsouth[pewdata$country=="Malaysia"] <- 1
pewdata$globalsouth[pewdata$country=="Mexico"] <- 1
pewdata$globalsouth[pewdata$country=="Netherlands"] <- 0
pewdata$globalsouth[pewdata$country=="Nigeria"] <- 1
pewdata$globalsouth[pewdata$country=="Philippines"] <- 1
pewdata$globalsouth[pewdata$country=="Poland"] <- 0
pewdata$globalsouth[pewdata$country=="Russia"] <- 0
pewdata$globalsouth[pewdata$country=="Senegal"] <- 1
pewdata$globalsouth[pewdata$country=="Singapore"] <- 0
pewdata$globalsouth[pewdata$country=="Slovakia"] <- 0
pewdata$globalsouth[pewdata$country=="South Africa"] <- 1
pewdata$globalsouth[pewdata$country=="South Korea"] <- 0
pewdata$globalsouth[pewdata$country=="Spain"] <- 0
pewdata$globalsouth[pewdata$country=="Sweden"] <- 0
pewdata$globalsouth[pewdata$country=="Tanzania"] <- 1
pewdata$globalsouth[pewdata$country=="Tunisia"] <- 1
pewdata$globalsouth[pewdata$country=="Turkey"] <- 0
pewdata$globalsouth[pewdata$country=="Uganda"] <- 1
pewdata$globalsouth[pewdata$country=="Ukraine"] <- 0
pewdata$globalsouth[pewdata$country=="United Kingdom"] <- 0
pewdata$globalsouth[pewdata$country=="United States"] <- 0
pewdata$globalsouth[pewdata$country=="Venezuela"] <- 1


pewdata$ally <- "" 
pewdata$ally[pewdata$country=="Argentina"] <- "MNNA"
pewdata$ally[pewdata$country=="Australia"] <- "MNNA"
pewdata$ally[pewdata$country=="Belgium"] <- "NATO"
pewdata$ally[pewdata$country=="Bolivia"] <- "Neutral"
pewdata$ally[pewdata$country=="Brazil"] <- "MNNA"
pewdata$ally[pewdata$country=="Bulgaria"] <- "NATO"
pewdata$ally[pewdata$country=="Canada"] <- "NATO"
pewdata$ally[pewdata$country=="Chile"] <- "Neutral"
pewdata$ally[pewdata$country=="Czech Republic"] <- "NATO"
pewdata$ally[pewdata$country=="Egypt"] <- "MNNA"
pewdata$ally[pewdata$country=="El Salvador"] <- "Neutral"
pewdata$ally[pewdata$country=="France"] <- "NATO"
pewdata$ally[pewdata$country=="Germany"] <- "NATO"
pewdata$ally[pewdata$country=="Ghana"] <- "Neutral"
pewdata$ally[pewdata$country=="Greece"] <- "NATO"
pewdata$ally[pewdata$country=="Hungary"] <- "NATO"
pewdata$ally[pewdata$country=="India"] <- "India- MDP"
pewdata$ally[pewdata$country=="Indonesia"] <- "Neutral"
pewdata$ally[pewdata$country=="Israel"] <- "MNNA"
pewdata$ally[pewdata$country=="Italy"] <- "NATO"
pewdata$ally[pewdata$country=="Japan"] <- "MNNA"
pewdata$ally[pewdata$country=="Jordan"] <- "MNNA"
pewdata$ally[pewdata$country=="Kenya"] <- "Neutral"
pewdata$ally[pewdata$country=="Lebanon"] <- "Neutral"
pewdata$ally[pewdata$country=="Lithuania"] <- "NATO"
pewdata$ally[pewdata$country=="Malaysia"] <- "Neutral"
pewdata$ally[pewdata$country=="Mexico"] <- "Neutral"
pewdata$ally[pewdata$country=="Netherlands"] <- "NATO"
pewdata$ally[pewdata$country=="Nigeria"] <- "Neutral"
pewdata$ally[pewdata$country=="Philippines"] <- "MNNA"
pewdata$ally[pewdata$country=="Poland"] <- "NATO"
pewdata$ally[pewdata$country=="Russia"] <- "Hostile"
pewdata$ally[pewdata$country=="Senegal"] <- "Neutral"
pewdata$ally[pewdata$country=="Singapore"] <- "Neutral"
pewdata$ally[pewdata$country=="Slovakia"] <- "NATO"
pewdata$ally[pewdata$country=="South Africa"] <- "Neutral"
pewdata$ally[pewdata$country=="South Korea"] <- "MNNA"
pewdata$ally[pewdata$country=="Spain"] <- "NATO"
pewdata$ally[pewdata$country=="Sweden"] <- "NATO"
pewdata$ally[pewdata$country=="Tanzania"] <- "Neutral"
pewdata$ally[pewdata$country=="Tunisia"] <- "MNNA"
pewdata$ally[pewdata$country=="Turkey"] <- "NATO"
pewdata$ally[pewdata$country=="Uganda"] <- "Neutral"
pewdata$ally[pewdata$country=="Ukraine"] <- "Ukraine"
pewdata$ally[pewdata$country=="United Kingdom"] <- "NATO"
pewdata$ally[pewdata$country=="United States"] <- "NATO"
pewdata$ally[pewdata$country=="Venezuela"] <- "Hostile"


# quick check that key DV is now scaled low to high
with(pewdata, table(country, demsat) )   # ok
# quick check that polarities are all in line with manual data audit
with(pewdata, cor(demsat, fav_us, use="complete.obs") ) # ok
with(pewdata, cor(demsat, fav_eu, use="complete.obs") ) # ok
with(pewdata, cor(fav_russia, fav_china, use="complete.obs") ) # ok
with(pewdata, cor(demsat, fav_un, use="complete.obs") ) # ok
with(pewdata, cor(demsat, fav_nato, use="complete.obs") ) # ok
with(pewdata, cor(fav_russia, confid_putin, use="complete.obs") ) # ok
with(pewdata, cor(fav_china, confid_xi, use="complete.obs") ) # ok
with(pewdata, cor(fav_us, confid_biden, use="complete.obs") ) # ok
with(pewdata, cor(fav_eu, confid_macron, use="complete.obs") ) # ok


# we put the whole original dataset into my.data
# this version of the file does not have 2023 (original submitted manuscript)
# we use the version with 2023 data for the new Appendix and some revised models. 

my.data <- read_dta("Pew2010-2022_VDEM.dta")

my.data$date <- as.Date(my.data$qdate_e)
table(my.data$date)

my.data$prewar <- NA 
my.data$prewar[my.data$date < as.Date("2022-02-25")] <- 1 
my.data$prewar[my.data$date > as.Date("2022-02-24")] <- 0 

my.data$postwar <- 1-my.data$prewar

my.data$sample <- 0 
my.data$sample[my.data$country=="Canada"] <- 1
my.data$sample[my.data$country=="France"] <- 1
my.data$sample[my.data$country=="Germany"] <- 1
my.data$sample[my.data$country=="Italy"] <- 1
my.data$sample[my.data$country=="Japan"] <- 1
my.data$sample[my.data$country=="United Kingdom"] <- 1
my.data$sample[my.data$country=="Sweden"] <- 1


my.data$globalsouth <- NA 

my.data$globalsouth[my.data$country=="Argentina"] <- 1
my.data$globalsouth[my.data$country=="Australia"] <- 0
# my.data$globalsouth[my.data$country=="Bangladesh"] <- 1
my.data$globalsouth[my.data$country=="Belgium"] <- 0
my.data$globalsouth[my.data$country=="Bolivia"] <- 1
my.data$globalsouth[my.data$country=="Brazil"] <- 1
my.data$globalsouth[my.data$country=="Britain"] <- 0
my.data$globalsouth[my.data$country=="Bulgaria"] <- 0
# my.data$globalsouth[my.data$country=="Burkina Faso"] <- 1
my.data$globalsouth[my.data$country=="Canada"] <- 0
my.data$globalsouth[my.data$country=="Chile"] <- 1
# my.data$globalsouth[my.data$country=="China"] <- 1
# my.data$globalsouth[my.data$country=="Colombia"] <- 1
my.data$globalsouth[my.data$country=="Czech Republic"] <- 0
# my.data$globalsouth[my.data$country=="Denmark"] <- 0
my.data$globalsouth[my.data$country=="Egypt"] <- 1
my.data$globalsouth[my.data$country=="El Salvador"] <- 1
# my.data$globalsouth[my.data$country=="Ethiopia"] <- 1
my.data$globalsouth[my.data$country=="France"] <- 0
my.data$globalsouth[my.data$country=="Germany"] <- 0
my.data$globalsouth[my.data$country=="Ghana"] <- 1
my.data$globalsouth[my.data$country=="Greece"] <- 0
my.data$globalsouth[my.data$country=="Hungary"] <- 0
my.data$globalsouth[my.data$country=="India"] <- 1
my.data$globalsouth[my.data$country=="Indonesia"] <- 1
my.data$globalsouth[my.data$country=="Israel"] <- 0
my.data$globalsouth[my.data$country=="Italy"] <- 0
my.data$globalsouth[my.data$country=="Japan"] <- 0
my.data$globalsouth[my.data$country=="Jordan"] <- 1
my.data$globalsouth[my.data$country=="Kenya"] <- 1
my.data$globalsouth[my.data$country=="Lebanon"] <- 1
my.data$globalsouth[my.data$country=="Lithuania"] <- 0
my.data$globalsouth[my.data$country=="Malaysia"] <- 1
my.data$globalsouth[my.data$country=="Mexico"] <- 1
my.data$globalsouth[my.data$country=="Netherlands"] <- 0
my.data$globalsouth[my.data$country=="Nigeria"] <- 1
my.data$globalsouth[my.data$country=="Philippines"] <- 1
my.data$globalsouth[my.data$country=="Poland"] <- 0
my.data$globalsouth[my.data$country=="Russia"] <- 0
my.data$globalsouth[my.data$country=="Senegal"] <- 1
my.data$globalsouth[my.data$country=="Singapore"] <- 0
my.data$globalsouth[my.data$country=="Slovakia"] <- 0
my.data$globalsouth[my.data$country=="South Africa"] <- 1
my.data$globalsouth[my.data$country=="South Korea"] <- 0
my.data$globalsouth[my.data$country=="Spain"] <- 0
my.data$globalsouth[my.data$country=="Sweden"] <- 0
my.data$globalsouth[my.data$country=="Tanzania"] <- 1
my.data$globalsouth[my.data$country=="Tunisia"] <- 1
my.data$globalsouth[my.data$country=="Turkey"] <- 0
my.data$globalsouth[my.data$country=="Uganda"] <- 1
my.data$globalsouth[my.data$country=="Ukraine"] <- 0
my.data$globalsouth[my.data$country=="United Kingdom"] <- 0
my.data$globalsouth[my.data$country=="United States"] <- 0
my.data$globalsouth[my.data$country=="Venezuela"] <- 1



my.data$ally <- "" 
my.data$ally[my.data$country=="Argentina"] <- "MNNA"
my.data$ally[my.data$country=="Australia"] <- "MNNA"
my.data$ally[my.data$country=="Belgium"] <- "NATO"
my.data$ally[my.data$country=="Bolivia"] <- "Neutral"
my.data$ally[my.data$country=="Brazil"] <- "MNNA"
my.data$ally[my.data$country=="Bulgaria"] <- "NATO"
my.data$ally[my.data$country=="Canada"] <- "NATO"
my.data$ally[my.data$country=="Chile"] <- "Neutral"
my.data$ally[my.data$country=="Czech Republic"] <- "NATO"
my.data$ally[my.data$country=="Egypt"] <- "MNNA"
my.data$ally[my.data$country=="El Salvador"] <- "Neutral"
my.data$ally[my.data$country=="France"] <- "NATO"
my.data$ally[my.data$country=="Germany"] <- "NATO"
my.data$ally[my.data$country=="Ghana"] <- "Neutral"
my.data$ally[my.data$country=="Greece"] <- "NATO"
my.data$ally[my.data$country=="Hungary"] <- "NATO"
my.data$ally[my.data$country=="India"] <- "India- MDP"
my.data$ally[my.data$country=="Indonesia"] <- "Neutral"
my.data$ally[my.data$country=="Israel"] <- "MNNA"
my.data$ally[my.data$country=="Italy"] <- "NATO"
my.data$ally[my.data$country=="Japan"] <- "MNNA"
my.data$ally[my.data$country=="Jordan"] <- "MNNA"
my.data$ally[my.data$country=="Kenya"] <- "Neutral"
my.data$ally[my.data$country=="Lebanon"] <- "Neutral"
my.data$ally[my.data$country=="Lithuania"] <- "NATO"
my.data$ally[my.data$country=="Malaysia"] <- "Neutral"
my.data$ally[my.data$country=="Mexico"] <- "Neutral"
my.data$ally[my.data$country=="Netherlands"] <- "NATO"
my.data$ally[my.data$country=="Nigeria"] <- "Neutral"
my.data$ally[my.data$country=="Philippines"] <- "MNNA"
my.data$ally[my.data$country=="Poland"] <- "NATO"
my.data$ally[my.data$country=="Russia"] <- "Hostile"
my.data$ally[my.data$country=="Senegal"] <- "Neutral"
my.data$ally[my.data$country=="Singapore"] <- "Neutral"
my.data$ally[my.data$country=="Slovakia"] <- "NATO"
my.data$ally[my.data$country=="South Africa"] <- "Neutral"
my.data$ally[my.data$country=="South Korea"] <- "MNNA"
my.data$ally[my.data$country=="Spain"] <- "NATO"
my.data$ally[my.data$country=="Sweden"] <- "NATO"
my.data$ally[my.data$country=="Tanzania"] <- "Neutral"
my.data$ally[my.data$country=="Tunisia"] <- "MNNA"
my.data$ally[my.data$country=="Turkey"] <- "NATO"
my.data$ally[my.data$country=="Uganda"] <- "Neutral"
my.data$ally[my.data$country=="Ukraine"] <- "Ukraine"
my.data$ally[my.data$country=="United Kingdom"] <- "NATO"
my.data$ally[my.data$country=="United States"] <- "NATO"
my.data$ally[my.data$country=="Venezuela"] <- "Hostile"

my.data$constant <- 0
my.data$constant[my.data$country=="Australia"] <- 1
my.data$constant[my.data$country=="Canada"] <- 1
my.data$constant[my.data$country=="France"] <- 1
my.data$constant[my.data$country=="Germany"] <- 1
my.data$constant[my.data$country=="Greece"] <- 1
my.data$constant[my.data$country=="Israel"] <- 1
my.data$constant[my.data$country=="Italy"] <- 1
my.data$constant[my.data$country=="Japan"] <- 1
my.data$constant[my.data$country=="Netherlands"] <- 1
my.data$constant[my.data$country=="Poland"] <- 1
my.data$constant[my.data$country=="South Korea"] <- 1
my.data$constant[my.data$country=="Spain"] <- 1
my.data$constant[my.data$country=="Sweden"] <- 1
my.data$constant[my.data$country=="United Kingdom"] <- 1

# these were scaled positive to negative, so for 4-point models we switch polarity
my.data$confid_putin <- -as.numeric(my.data$confid_putin)
my.data$confid_xi <- -as.numeric(my.data$confid_xi)
my.data$confid_uspres <- rowMeans( cbind(as.numeric(my.data$confid_trump), as.numeric(my.data$confid_biden), as.numeric(my.data$confid_obama)), na.rm=T )
my.data$confid_uspres <- -as.numeric(my.data$confid_uspres)
my.data$confid_biden <-  -(as.numeric(my.data$confid_biden))
my.data$confid_macron <-  -(as.numeric(my.data$confid_macron))
my.data$confid_scholz <-  -(as.numeric(my.data$confid_scholz))
my.data$confid_baseline <- - as.numeric(my.data$confid_baseline)


# code some dummy variables for 4-point scale items

my.data$confid_biden_01 <- NA 
my.data$confid_biden_01[my.data$confid_biden==0] <- 0 
my.data$confid_biden_01[my.data$confid_biden==1] <- 0 
my.data$confid_biden_01[my.data$confid_biden==2] <- 1 
my.data$confid_biden_01[my.data$confid_biden==3] <- 1 

my.data$confid_xi_01 <- NA 
my.data$confid_xi_01[my.data$confid_xi==0] <- 0 
my.data$confid_xi_01[my.data$confid_xi==1] <- 0 
my.data$confid_xi_01[my.data$confid_xi==2] <- 1 
my.data$confid_xi_01[my.data$confid_xi==3] <- 1 

my.data$confid_putin_01 <- NA 
my.data$confid_putin_01[my.data$confid_putin==0] <- 0 
my.data$confid_putin_01[my.data$confid_putin==1] <- 0 
my.data$confid_putin_01[my.data$confid_putin==2] <- 1 
my.data$confid_putin_01[my.data$confid_putin==3] <- 1 

my.data$fav_us_01 <- NA 
my.data$fav_us_01[my.data$fav_us==1] <- 0 
my.data$fav_us_01[my.data$fav_us==2] <- 0 
my.data$fav_us_01[my.data$fav_us==3] <- 1 
my.data$fav_us_01[my.data$fav_us==4] <- 1 


my.data$fav_nato <- -as.numeric(my.data$fav_nato)
my.data$fav_un <- -as.numeric(my.data$fav_un)


# quick check that now scaled low to high 
with(my.data, table(country, demsat) )
# quick check that polarities are all in line with manual data audit
with(my.data, cor(demsat, fav_us, use="complete.obs") ) # ok
with(my.data, cor(demsat, fav_eu, use="complete.obs") ) # ok
with(my.data, cor(fav_russia, fav_china, use="complete.obs") ) # ok
with(my.data, cor(demsat, fav_un, use="complete.obs") ) # ok
with(my.data, cor(demsat, fav_nato, use="complete.obs") ) # ok
with(my.data, cor(fav_russia, confid_putin, use="complete.obs") ) # ok
with(my.data, cor(fav_china, confid_xi, use="complete.obs") ) # ok
with(my.data, cor(fav_us, confid_biden, use="complete.obs") ) # ok
with(my.data, cor(fav_eu, confid_macron, use="complete.obs") ) # ok

# confirm correct recoding for dummies 
with(my.data, cor(confid_putin, confid_putin_01, use="complete.obs") ) # ok
with(my.data, cor(confid_xi, confid_xi_01, use="complete.obs") ) # ok
with(my.data, cor(confid_biden, confid_biden_01, use="complete.obs") ) # ok

# check prowestern variable is still correctly scaled since generated in Stata
with(my.data, cor(prowestern, fav_us, use="complete.obs") ) # ok
with(my.data, cor(prowestern, fav_russia, use="complete.obs") ) # ok



## We use this data frame just for the 2022 subset. 
## It is the same as my.data for 2022 only, with some recoding done already in Stata.

data2022 <-  read_dta("Pew2022recoded_experiment.dta")

data2022$country <- as.character(as.factor(data2022$country))

# cbeck coding of variables and completeness
table(data2022$country)
table(data2022$date)

# pre or post dummy variable 
data2022$prewar <- NA 
data2022$prewar[data2022$date < as.Date("2022-02-25")] <- 1 
data2022$prewar[data2022$date > as.Date("2022-02-24")] <- 0 

data2022$postwar <- 1-data2022$prewar

# we are still coding Sweden as part of the sample here, but exclude it later,
# due to low prewar sample (polling began only day before invasion)

data2022$sample <- 0 
data2022$sample[data2022$country=="Canada"] <- 1 ## canada
data2022$sample[data2022$country=="France"] <- 1 ## france
data2022$sample[data2022$country=="Germany"] <- 1 ## germany
data2022$sample[data2022$country=="Italy"] <- 1 ## italy
data2022$sample[data2022$country=="Japan"] <- 1 ## japan
data2022$sample[data2022$country=="UK"] <- 1 ## UK
data2022$sample[data2022$country=="Sweden"] <- 1 ## Sweden

# check codings
table(data2022$sample)
table(subset(data2022, date < as.Date("2022-02-25"))$sample)
table(subset(data2022, date < as.Date("2022-02-25"))$country)
table(subset(data2022, date < as.Date("2022-02-24"))$country)
table(data2022$citizen_vote)

# just in case residual NA codings here
#data2022$citizen_vote[data2022$citizen_vote==8] <- NA
#data2022$citizen_vote[data2022$citizen_vote==9] <- NA

# adding coding for this variable
data2022$vote_important <- NA 
data2022$vote_important[data2022$citizen_vote==1] <- 0 
data2022$vote_important[data2022$citizen_vote==2] <- 0 
data2022$vote_important[data2022$citizen_vote==3] <- 1 
data2022$vote_important[data2022$citizen_vote==4] <- 1 

# in case missing 
# data2022$age_square <- data2022$age^2
# data2022$female <- data2022$sex

# again these ones were missed in initial manual polarity recoding 
# as were not used in initial article draft/models 
data2022$fav_un <- 5-data2022$fav_un
data2022$fav_nato <- 5-data2022$fav_nato
data2022$confid_putin <- 5-data2022$confid_putin
data2022$confid_xi <- 5-data2022$confid_xi
data2022$confid_biden <- 5-data2022$confid_biden
data2022$confid_macron <- 5-data2022$confid_macron
data2022$confid_scholz <- 5-data2022$confid_scholz


# quick check that now scaled low to high 
with(data2022, table(country, demsat) ) # ok
# quick check that polarities are all in line with manual data audit
with(data2022, cor(demsat, fav_us, use="complete.obs") ) # ok
with(data2022, cor(demsat, fav_eu, use="complete.obs") ) # ok
with(data2022, cor(fav_russia, fav_china, use="complete.obs") ) # ok
with(data2022, cor(demsat, fav_un, use="complete.obs") ) # ok
with(data2022, cor(demsat, fav_nato, use="complete.obs") ) # ok
with(data2022, cor(fav_russia, confid_putin, use="complete.obs") ) # ok
with(data2022, cor(fav_china, confid_xi, use="complete.obs") ) # ok
with(data2022, cor(fav_us, confid_biden, use="complete.obs") ) # ok
with(data2022, cor(fav_eu, confid_macron, use="complete.obs") ) # ok


## MODELS AND FIGURES

## Figure 1 - cross-country favorability comparison chart. 

## This data comes from a separate project, and is imported specifically for this figure. 

ddata <- read.csv("ddata.csv")
global_series <- read.csv("global_series.csv")


ld2 <- 
  ggplot( subset(ddata, libdemoc==1 & year==2022), aes(x=china_combined*100,y=russia_combined_series)) + 
  geom_point(data=subset(ddata, libdemoc==1 & year==2022), aes(x=china_combined*100,y=russia_combined_series),
             size=sqrt(subset(ddata, libdemoc==1 & year==2022)$pop)/1000 , alpha=0.75) + 
  geom_text_repel(data=subset(global_series, year==2022), aes(x=china*100, y=russia, label="Global\nAverage"), nudge_x=30, nudge_y=0, alpha=0.25) + 
  annotate(geom="point", x=subset(global_series, year==2022)$china*100, y=subset(global_series, year==2022)$russia, size=6, color="black"  ) + 
  annotate(geom="point", x=subset(global_series, year==2022)$china*100, y=subset(global_series, year==2022)$russia, size=5, color="white"  ) + 
  theme_article(base_size=16) + 
  theme(panel.grid.major.y = element_blank()) + 
  theme(panel.grid.minor.y = element_blank()) + 
  theme(panel.border = element_blank()) +
  theme(plot.background = element_blank()) + 
  coord_cartesian(xlim=c(3,100), ylim=c(3,100)) + 
  xlab("Favorable towards China (%)")+ 
  ylab("Favorable towards Russia (%)") + 
  scale_y_continuous(labels =function(x) paste0(x,"%")) + 
  scale_x_continuous(labels =function(x) paste0(x,"%")) +
  theme(plot.title = element_text(hjust = 0.5) ) + 
  ggtitle("2022")


ld1 <- 
  ggplot( subset(ddata, libdemoc==1 & year==2012), aes(x=china_combined*100,y=russia_combined_series)) + 
  geom_point(data=subset(ddata, libdemoc==1 & year==2012), aes(x=china_combined*100,y=russia_combined_series),
             size=sqrt(subset(ddata, libdemoc==1 & year==2012)$pop)/1000 , alpha=0.75) + 
  geom_text_repel(data=subset(global_series, year==2012), aes(x=china*100, y=russia, label="Global\nAverage"), nudge_x=30, nudge_y=0, alpha=0.25) + 
  annotate(geom="point", x=subset(global_series, year==2012)$china*100, y=subset(global_series, year==2012)$russia, size=6, color="black"  ) + 
  annotate(geom="point", x=subset(global_series, year==2012)$china*100, y=subset(global_series, year==2012)$russia, size=5, color="white"  ) + 
  theme_article(base_size=16) + 
  theme(panel.grid.major.y = element_blank()) + 
  theme(panel.grid.minor.y = element_blank())+ 
  theme(panel.border = element_blank()) +
  theme(plot.background = element_blank()) + 
  coord_cartesian(xlim=c(3,100), ylim=c(3,100))+ 
  xlab("Favorable towards China (%)")+ 
  ylab("Favorable towards Russia (%)") + 
  scale_y_continuous(labels =function(x) paste0(x,"%")) + 
  scale_x_continuous(labels =function(x) paste0(x,"%")) + 
  theme(plot.title = element_text(hjust = 0.5) ) + 
  ggtitle("2012")


library(cowplot)

libdemoc_ruschina_scatter <- plot_grid(ld1, ld2)

# save.plot("libdemoc_ruschina_scatter.png", width=10, height=5)

save.tiff(libdemoc_ruschina_scatter,
          "libdemoc_ruschina_scatter.tiff",
          w=10, h=5)


# models for pre and post invasion coefficients.
# we can use my.data for this, as the postwar variable automatically subsets to 2022 only. 

# valid entries for !is.na(postwar) and sample==1 are only the 2022 sample countries. 
table(subset(my.data, !is.na(postwar))$country)
table(subset(my.data, sample==1)$country)

# in the order of Table 1


# Favor Russia 
exp_6a <-  lm(data=subset(my.data, sample==1), fav_russia ~ postwar + factor(country) , weights=weight )
exp_6b <-  lm(data=subset(my.data, sample==1), fav_russia ~ postwar + age + age_square + female + factor(country) + education + income , weights=weight )
exp_6c <-  lm(data=subset(my.data, sample==1), fav_russia ~ postwar + age + age_square + female + factor(country) + education + income + fav_baseline , weights=weight )

# Favor China
exp_7a <-  lm(data=subset(my.data, sample==1), fav_china ~ postwar + factor(country) , weights=weight )
exp_7b <-  lm(data=subset(my.data, sample==1), fav_china ~ postwar + age + age_square + female + factor(country) + education + income , weights=weight )
exp_7c <-  lm(data=subset(my.data, sample==1), fav_china ~ postwar + age + age_square + female + factor(country) + education + income + fav_baseline, weights=weight )

# Favor US
exp_8a <-  lm(data=subset(my.data, sample==1), fav_us ~ postwar + factor(country), weights=weight )
exp_8b <-  lm(data=subset(my.data, sample==1), fav_us ~ postwar + age + age_square + female + factor(country) + education + income  , weights=weight )
exp_8c <-  lm(data=subset(my.data, sample==1), fav_us ~ postwar + age + age_square + female + factor(country) + education + income + fav_baseline  , weights=weight )


# Putin Approval
exp_3a <- lm(data=subset(my.data, sample==1), confid_putin ~ postwar + factor(country) , weights=weight )
exp_3b <- lm(data=subset(my.data, sample==1), confid_putin ~ postwar + age + age_square + female + factor(country) + education + income , weights=weight )
exp_3c <- lm(data=subset(my.data, sample==1), confid_putin ~ postwar + age + age_square + female + factor(country) + education + income + confid_baseline, weights=weight )

# Xi Approval
exp_2a <- lm(data=subset(my.data, sample==1), confid_xi ~ postwar + factor(country) , weights=weight )
exp_2b <- lm(data=subset(my.data, sample==1), confid_xi ~ postwar + age + age_square + female + factor(country) + education + income, weights=weight )
exp_2c <- lm(data=subset(my.data, sample==1), confid_xi ~ postwar + age + age_square + female + factor(country) + education + income  + confid_baseline, weights=weight )

# Biden Approval
exp_1a <- lm(data=subset(my.data, sample==1), confid_biden ~ postwar +  factor(country)  , weights=weight )
exp_1b <- lm(data=subset(my.data, sample==1), confid_biden ~ postwar + age + age_square + female + factor(country) + education  + income , weights=weight )
exp_1c <- lm(data=subset(my.data, sample==1), confid_biden ~ postwar + age + age_square + female + factor(country) + education  + income + confid_baseline , weights=weight )

# Macron Approval
exp_4a <- lm(data=subset(my.data, sample==1), confid_macron ~ postwar + factor(country), weights=weight )
exp_4b <- lm(data=subset(my.data, sample==1), confid_macron ~ postwar + age + age_square + female + factor(country) + education + income, weights=weight )
exp_4c <- lm(data=subset(my.data, sample==1), confid_macron ~ postwar + age + age_square + female + factor(country) + education + income + confid_baseline, weights=weight )

# Scholz Approval
exp_5a <- lm(data=subset(my.data, sample==1), confid_scholz ~ postwar + factor(country), weights=weight )
exp_5b <- lm(data=subset(my.data, sample==1), confid_scholz ~ postwar + age + age_square + female + factor(country) + education + income , weights=weight )
exp_5c <- lm(data=subset(my.data, sample==1), confid_scholz ~ postwar + age + age_square +  female + factor(country) + education + income + confid_baseline, weights=weight )

## China military a problem: again polarity check. 
my.data$china_milpower <- 5 - as.numeric(my.data$china_milpower)
# check polarity is correct: should now be negative
with(my.data, cor(fav_china, china_milpower, use="complete.obs") ) # ok
with(my.data, cor(confid_xi, china_milpower, use="complete.obs") ) # ok


exp_12a <-  lm(data=subset(my.data, sample==1), china_milpower ~ postwar + factor(country), weights=weight )
exp_12b <-  lm(data=subset(my.data, sample==1), china_milpower ~ postwar + age + age_square + female + factor(country) + education + income  , weights=weight )

## US reliable partner variable: needs polarity check.
my.data$reliable_us <- 5 - as.numeric(my.data$reliable_us)
# check polarity is correct: should now be positive
with(my.data, cor(fav_us, reliable_us, use="complete.obs") ) # ok
with(my.data, cor(confid_biden, reliable_us, use="complete.obs") ) # ok


exp_11a <-  lm(data=subset(my.data, sample==1), reliable_us ~ postwar + factor(country), weights=weight )
exp_11b <-  lm(data=subset(my.data, sample==1), reliable_us ~ postwar + age + age_square + female + factor(country) + education + income  , weights=weight )

# voting important to be good member of society. 
# adding this to the model required going back and getting a version of the data with this extra variable.
exp_14a <-  lm(data=subset(data2022, sample==1), vote_important ~ postwar + factor(country), weights=weight)
exp_14b <-  lm(data=subset(data2022, sample==1), vote_important ~ postwar + female + age + age_square + income + education + factor(country), weights=weight)


## output table for models.
## we only took the postwar row for each for parsimony.
## these were combined direct in TexMaker.

stargazer(exp_6a, exp_6b, exp_6c,
          exp_7a, exp_7b, exp_7c)
stargazer(exp_8a, exp_8b, exp_8c,
          exp_3a, exp_3b, exp_3c)
stargazer(exp_2a, exp_2b, exp_2c,
          exp_1a, exp_1b, exp_1c) 
stargazer(exp_4a, exp_4b, exp_4c,
          exp_5a, exp_5b, exp_5c)
stargazer(exp_12a, exp_12b, 
          exp_11a, exp_11b, 
          exp_14a, exp_14b)  


## Figure 2: Effect of democratic legitimacy, pre and post. 
## Sweden manually excluded here. 
## Doesn't make huge difference to estimates but best to be safe, 
## and avoid objections over fact that Sweden polling began so close to discontinuity date.


ds1a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  confid_xi ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds1b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  confid_xi ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds2a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  confid_biden ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds2b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  confid_biden ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds3a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  confid_putin ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds3b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  confid_putin ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds4a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  confid_macron ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds4b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  confid_macron ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds5a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  confid_scholz ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds5b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  confid_scholz ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds6a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  fav_eu ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds6b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  fav_eu ~ demsat + age + age_square + female + factor(country) + education + income  
           , weights=weight )

ds7a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  fav_china ~ demsat + age + age_square + female + factor(country) + education + income 
           , weights=weight )

ds7b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  fav_china ~ demsat + age + age_square + female + factor(country) + education + income  
           , weights=weight )

ds8a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==0),  fav_russia ~ demsat + age + age_square + female + factor(country) + education + income  
           , weights=weight )

ds8b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & postwar==1),  fav_russia ~ demsat + age + age_square + female + factor(country) + education + income  
           , weights=weight )

ds9a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & date<as.Date("2022-02-25")),  fav_us ~ demsat + age + age_square + female + factor(country) + education + income
           , weights=weight )

ds9b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & date>as.Date("2022-02-24")),  fav_us ~ demsat + age + age_square + female + factor(country) + education + income  
           , weights=weight )

ds10a <- lm(data=subset(my.data, sample==1 & country!="Sweden" & date<as.Date("2022-02-25")),  fav_nato ~ demsat + age + age_square + female + factor(country) + education + income 
            , weights=weight )

ds10b <- lm(data=subset(my.data, sample==1 & country!="Sweden" & date>as.Date("2022-02-24")),  fav_nato ~ demsat + age + age_square + female + factor(country) + education + income
            , weights=weight )


# make a nice data frame table for these. 

results <- 
  as.data.frame (matrix(data=c(NA,NA,NA,NA,NA), nrow=5, ncol=5) )

colnames(results) <- c("var", "pre", "post", "pre_sd", "post_sd")

results$var[1] <- "Favorable to European Union"
results$var[2] <- "Favorable to NATO"
results$var[3] <- "Favorable to United States"
results$var[4] <- "Favorable to China"
results$var[5] <- "Favorable to Russia"

results$pre[1] <- summary(ds6a)$coefficients["demsat","Estimate"]
results$post[1] <- summary(ds6b)$coefficients["demsat","Estimate"]
results$pre_sd[1] <- summary(ds6a)$coefficients["demsat","Std. Error"]
results$post_sd[1] <- summary(ds6b)$coefficients["demsat","Std. Error"]

results$pre[2] <- summary(ds10a)$coefficients["demsat","Estimate"]
results$post[2] <- summary(ds10b)$coefficients["demsat","Estimate"]
results$pre_sd[2] <- summary(ds10a)$coefficients["demsat","Std. Error"]
results$post_sd[2] <- summary(ds10b)$coefficients["demsat","Std. Error"]

results$pre[3] <- summary(ds9a)$coefficients["demsat","Estimate"]
results$post[3] <- summary(ds9b)$coefficients["demsat","Estimate"]
results$pre_sd[3] <- summary(ds9a)$coefficients["demsat","Std. Error"]
results$post_sd[3] <- summary(ds9b)$coefficients["demsat","Std. Error"]

results$pre[4] <- summary(ds7a)$coefficients["demsat","Estimate"]
results$post[4] <- summary(ds7b)$coefficients["demsat","Estimate"]
results$pre_sd[4] <- summary(ds7a)$coefficients["demsat","Std. Error"]
results$post_sd[4] <- summary(ds7b)$coefficients["demsat","Std. Error"]

results$pre[5] <- summary(ds8a)$coefficients["demsat","Estimate"]
results$post[5] <- summary(ds8b)$coefficients["demsat","Estimate"]
results$pre_sd[5] <- summary(ds8a)$coefficients["demsat","Std. Error"]
results$post_sd[5] <- summary(ds8b)$coefficients["demsat","Std. Error"]

results$var <- factor(results$var, levels=c("Favorable to Russia","Favorable to China",
                                            "Favorable to United States",
                                            "Favorable to NATO","Favorable to European Union"))


demsat_shift <- 
ggplot(data=results, aes(x=var, y=pre) ) + 
  geom_point(size=6, color="darkgrey") + 
  geom_segment(data=results, aes(x=var, xend=var, y=post-1.645*post_sd, yend=post+1.645*post_sd), color="black") +
  geom_segment(data=subset(results, post < pre-0.035),
               aes(x=var, xend=var, y=pre, yend=post+0.03), size=1.5, 
               color="darkgrey", arrow = arrow(type = "closed",length = unit(5, "points")) ) +
  geom_segment(data=subset(results, post > pre+0.025), size=1.5,
               aes(x=var, xend=var, y=pre, yend=post-0.025), color="darkgrey", arrow = arrow(type = "closed",length = unit(5, "points")) ) +
  geom_segment(data=subset(results, var=="Favorable to United States"), size=1.5,
               aes(x=var, xend=var, y=pre, yend=post-0.01), color="darkgrey", arrow = arrow(type = "closed",length = unit(5, "points")) ) +
  geom_text(data=subset(results, post < pre & var=="Favorable to China"), label="Pre-Invasion Sample", 
            color="#888888", hjust=0, size=3.4, aes(x=var, y=pre+0.0185)) +
  geom_text(data=subset(results, post < pre & var=="Favorable to China"), label="Post", 
            color="#444444", hjust=1, size=3.4, aes(x=var, y=post-0.025)) +
  geom_text(data=subset(results, post > pre & var=="Favorable to NATO"), label="Pre-Invasion Sample", 
            color="#888888", hjust=1, size=3.4, aes(x=var, y=pre-0.0185)) +
  geom_text(data=subset(results, post > pre & var=="Favorable to NATO"), label="Post", 
            color="#444444", hjust=0, size=3.4, aes(x=var, y=post+0.025)) +
  geom_point(size=2, color="white") + 
  geom_point(data=results, aes(x=var, y=post), size=6, color="black" ) + 
  geom_point(data=results, aes(x=var, y=post), size=2, color="white" ) + 
  theme_article() + 
  theme(panel.grid.major.y=element_blank()) +
  coord_flip(ylim=c(-0.15, 0.45)) + 
  xlab("") + ylab("Satisfaction with Democracy Coefficient") 


#save.plot("demsat_shift.png", width=8, height=2.5)

# for final version we want a tiff version 

save.tiff(demsat_shift, "demsat_shift.tiff", w=8, h=2.5)


## ORDERED LOGIT MODELS 

## For the final manuscript, we re-estimated adding the 2023 data. 

## These were done in Stata; see separate .do file. 



## MULTILEVEL MODELS 

## In the final paper, we are using pewdata for this:
## This is because we revised to include 2023 results as well in the models. 


## Country favorability over time slopes (Figure 3)

## constant country sample. 
## here we want to re-estimate MLM time series to ensure only countries 
## which are included on a constant(ish) basis, to prevent sample changes
## from biasing our coefficients. 

pewdata$constant <- 0
pewdata$constant[pewdata$country=="Australia"] <- 1
pewdata$constant[pewdata$country=="Canada"] <- 1
pewdata$constant[pewdata$country=="France"] <- 1
pewdata$constant[pewdata$country=="Germany"] <- 1
pewdata$constant[pewdata$country=="Greece"] <- 1
pewdata$constant[pewdata$country=="Israel"] <- 1
pewdata$constant[pewdata$country=="Italy"] <- 1
pewdata$constant[pewdata$country=="Japan"] <- 1
pewdata$constant[pewdata$country=="Netherlands"] <- 1
pewdata$constant[pewdata$country=="Poland"] <- 1
pewdata$constant[pewdata$country=="South Korea"] <- 1
pewdata$constant[pewdata$country=="Spain"] <- 1
pewdata$constant[pewdata$country=="Sweden"] <- 1
pewdata$constant[pewdata$country=="United Kingdom"] <- 1

## There was also a robustness check with an "expanded" constant sample.
## "expanded" meant countries that dropped out for a few years, but then came back. 

pewdata$constant_plus <- pewdata$constant
pewdata$constant_plus[pewdata$country=="Hungary"] <- 1
pewdata$constant_plus[pewdata$country=="Argentina"] <- 1
pewdata$constant_plus[pewdata$country=="Brazil"] <- 1
pewdata$constant_plus[pewdata$country=="India"] <- 1
pewdata$constant_plus[pewdata$country=="Indonesia"] <- 1
pewdata$constant_plus[pewdata$country=="Kenya"] <- 1
pewdata$constant_plus[pewdata$country=="Mexico"] <- 1
pewdata$constant_plus[pewdata$country=="Nigeria"] <- 1
pewdata$constant_plus[pewdata$country=="South Africa"] <- 1


# estimate the basic MLM, this time just with constant sample. 
# We estimate each then get slopes and coefs in a new frame and the 90% cis.

ml_12_1 <- lmer(data= subset(pewdata, constant==1), fav_china ~   age + age_square + female  + factor(country) + education + income + econ_sit + fav_baseline + 
                  ( 1 + demsat | year ), weights=weight )

grue1 <- ranef(ml_12_1)$year
grue1 <- adf(grue1)
grue1$year <- as.numeric( row.names(grue1) )

grue1$se.demsat <- se.ranef(ml_12_1)$year[,"demsat"]
grue1$demsat.lower <- grue1$demsat - 1.645*grue1$se.demsat
grue1$demsat.upper <- grue1$demsat + 1.645*grue1$se.demsat

ml_12_2 <- lmer(data= subset(pewdata, constant==1), fav_russia  ~ age + age_square + female  + factor(country) +  education + income + fav_baseline + econ_sit + 
                  ( 1 + demsat | year ), weights=weight )

grue2 <- ranef(ml_12_2)$year
grue2 <- adf(grue2)
grue2$year <- as.numeric( row.names(grue2) )

grue2$se.demsat <- se.ranef(ml_12_2)$year[,"demsat"]
grue2$demsat.lower <- grue2$demsat - 1.645*grue2$se.demsat
grue2$demsat.upper <- grue2$demsat + 1.645*grue2$se.demsat

ml_12_3 <- lmer(data= subset(pewdata, constant==1), fav_us  ~ age + age_square + female  + factor(country) + education + income +  econ_sit + fav_baseline + 
                  ( 1 + demsat | year ), weights=weight )

grue3 <- ranef(ml_12_3)$year
grue3 <- adf(grue3)
grue3$year <- as.numeric( row.names(grue3) )

grue3$se.demsat <- se.ranef(ml_12_3)$year[,"demsat"]
grue3$demsat.lower <- grue3$demsat - 1.645*grue3$se.demsat
grue3$demsat.upper <- grue3$demsat + 1.645*grue3$se.demsat

ml_12_4 <- lmer(data= subset(pewdata, constant==1), fav_eu  ~ age + age_square + female  + factor(country) + education + income +  econ_sit + fav_baseline + 
                  ( 1 + demsat | year ), weights=weight )

grue4 <- ranef(ml_12_4)$year
grue4 <- adf(grue4)
grue4$year <- as.numeric( row.names(grue4) )

grue4$se.demsat <- se.ranef(ml_12_4)$year[,"demsat"]
grue4$demsat.lower <- grue4$demsat - 1.645*grue4$se.demsat
grue4$demsat.upper <- grue4$demsat + 1.645*grue4$se.demsat

# Note: 2017 omitted due to country sample bias
MLM_slopes_1 <- 
ggplot(subset(grue1, year>2013 & year!=2018), aes(x=year, y=demsat)) +
  annotate(geom="rect", xmin=2021, xmax=2022, ymin=-Inf, ymax=Inf, fill="grey", alpha=0.35) +
  geom_vline(xintercept=2022, color="grey", linetype="dashed", size=0.35)+
  geom_vline(xintercept=2021, color="grey", linetype="dashed", size=0.35)+
  geom_ribbon(aes(ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_ribbon(data=subset(grue2,year>2013& year!=2017& year!=2018), aes(x=year, ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_ribbon(data=subset(grue3,year>2013& year!=2017& year!=2018), aes(x=year, ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_ribbon(data=subset(grue4, year>2013& year!=2017& year!=2018), aes(x=year, ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_segment(data= subset(grue1, year==2023 | year==2014), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data= subset(grue2, year==2023 | year==2014), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data= subset(grue3, year==2023 | year==2014), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data= subset(grue4, year==2023 | year==2014), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_line(color="darkgrey", size=2) + 
  geom_line(data=subset(grue2, year>2013& year!=2017& year!=2018), aes(year, demsat), size=2) + 
  geom_line(data=subset(grue3, year>2013& year!=2017& year!=2018), aes(year, demsat  ), size=2, color="#444444") +
  geom_line(data=subset(grue4, year>2013& year!=2017& year!=2018), aes(year, demsat  ), size=2)  +
  geom_line(data=subset(grue4, year>2013& year!=2017& year!=2018), aes(year, demsat  ), size=1, color="white")  +
  geom_point(data=subset(grue1, year==2023 | year==2014), aes(year,demsat)) + 
  geom_point(data=subset(grue2, year==2023 | year==2014), aes(year,demsat)) + 
  geom_point(data=subset(grue3, year==2023 | year==2014), aes(year,demsat)) + 
  geom_point(data=subset(grue4, year==2023 | year==2014), aes(year,demsat)) + 
  xlab("") + 
  ylab("Satisfaction with Democracy Estimated Effect") + 
  theme_article(base_size=15) +  
  scale_x_continuous( breaks=c(2010, 2012,2014,2016,2018,2020,2022), limits=c(2012,2024)) + 
  scale_y_continuous( breaks=c(-0.2,-0.1,0,0.1,0.2,0.3,0.4), labels=c("-0.20","-0.10","0.00","+0.10","+0.20","+0.30","+0.40")) + 
  coord_cartesian(xlim=c(2013.5,2025), ylim=c(-0.05, 0.31)) + 
  annotate(geom="text", label="Pre/Post\n Ukraine \nInvasion", hjust=1, x=2021.9, y= 0.305, size=3.3 ) +
  annotate(geom="text", label="Favor China", hjust=0, x=2023.2, y= tail(grue1$demsat,1) ) +
  annotate(geom="text", label="Favor Russia", hjust=0, x=2023.2, y= tail(grue2$demsat,1)  ) +
  annotate(geom="text", label="Favor the US", hjust=0, x=2023.2, y= tail(grue3$demsat,1)  ) +
  annotate(geom="text", label="Favor the EU", hjust=0, x=2023.2, y= tail(grue4$demsat,1)  ) 


# save.plot("MLM_slopes_demsat_bycountry_2014_2023_ccsample.png", width=9)

save.tiff(MLM_slopes_1, "MLM_slopes_demsat_bycountry_2014_2023_ccsample.tiff",
          w=9, h=6)


## Figure 4, MLM Slopes by Country 

ml_31 <- lmer(data=pewdata, prowestern ~ age + age_square + female + factor(country) + factor(year) + education + income + religion_import +  
                ( 1 + demsat | country ) , weights=weight )

soo <- ranef(ml_31)$country
soo <- adf(soo)
soo

soo$country <- row.names(soo)

vdem <- read.dta13("libdem.dta")

vdem <- subset(vdem, year==2019)
vdem$country[vdem$country=="United States of America"] <- "United States"
soo <- merge(soo, vdem, by=c("country"), all.x=T, all.y=F)

## ggrepel not working well for the .tiff export, so need to add a manual offset variable. 

soo$nudge.x <- NA
soo$nudge.x[soo$demsat>0.1] <- 0.02
soo$nudge.x[soo$demsat<0.1] <- -0.02
soo$nudge.y <- NA
soo$nudge.y[soo$v2x_libdem>0.36] <- 0.025
soo$nudge.y[soo$v2x_libdem<0.36] <- -0.025

libdem_slopes <- 
ggplot(soo, aes(demsat, v2x_libdem)) + 
  geom_smooth(method="lm", size=2, color="black", fill="#DDDDDD") + 
  geom_point(size=3) + 
  geom_point(size=1, color="white") + 
  geom_text_repel(data=soo, aes(label=country), hjust = "outward", nudge_x  = soo$nudge.x,
                  nudge_y=soo$nudge.y, alpha=0.65) + 
  theme_article() + 
  coord_cartesian(ylim=c(0,0.9), xlim=c(-1,1) ) +
  xlab("Country Random Slope, Democratic Satisfaction on Pro-Western Attitudes") + 
  ylab("Liberal Democracy Index")

# save.plot("country_slopes_vdem.png", width=7, height=7)

save.tiff(libdem_slopes, "country_slopes_vdem.tiff",
          w=7, h=7)


## Figure 5, by Alliance Status 

ml_33 <- lmer(data=pewdata, prowestern ~ age + age_square + female + factor(year) + education + income + religion_import +  
                ( 1 + demsat | ally ), weights=weight  )

summary(ml_33)

boo <- ranef(ml_33)$ally
boo <- adf(boo)
boo

boo$se.demsat <- se.ranef(ml_33)$ally[,"demsat"]
boo$demsat.lower <- boo$demsat - 1.645*boo$se.demsat
boo$demsat.upper <- boo$demsat + 1.645*boo$se.demsat

boo$name <-  row.names(boo) 
boo$name <- c("Adverse States", "India (Major Defense Partner)", "Major Non-NATO Allies", "NATO Allies", "Neutral States", "Ukraine")
boo$name <- as.factor(boo$name)

boo <- stort(boo, -boo$demsat)
boo$name <- factor(boo$name,levels= c("Adverse States", "Neutral States", "India (Major Defense Partner)", "Major Non-NATO Allies", "NATO Allies",  "Ukraine") )

alliance_effects <- 
ggplot(data=boo, aes(name, demsat)) + 
  geom_segment(aes(x=name,xend=name,y=-Inf,yend=Inf), linetype="dashed", color="white") + 
  geom_hline(yintercept=0) + 
  geom_segment(aes(x=name,xend=name,y=demsat.lower,yend=demsat.upper)) + 
  geom_point(size=3) + 
  geom_point(size=2, color="white") + 
  coord_flip() + 
  theme_article(base_size=13.25) + 
  xlab("") + ylab("Estimated Effect - Satisfaction with Democracy, Pro-Western Orientation")

# save.plot("Alliance_Effects.png", width=9.25, height=4)

save.tiff(alliance_effects, "Alliance_Effects.tiff",
          w=9.25, h=4)



## APPENDIX 
## additional tables and figures for the Appendix. 
## mainly descriptive tables, some T-tests and other minor requested extras.  

## plot for the new Appendix. 

plot.data <- as.data.frame(cbind(c("Favor Russia","Favor China","Favor the US",
                                   "Confidence in Biden","Confidence in Macron",
                                   "Confidence in Scholz","Confidence in Xi","Confidence in Putin"), 
                                 c(NA,NA,NA,NA,NA,NA,NA,NA),
                                 c(NA,NA,NA,NA,NA,NA,NA,NA),
                                 c(NA,NA,NA,NA,NA,NA,NA,NA) ))

names(plot.data) <- c("Variable", "Coefficient", "Lower", "Upper")

plot.data$Coefficient <- as.numeric(plot.data$Coefficient)


plot.data$Coefficient[1] <- as.numeric(summary(exp_6a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[1] <- plot.data$Coefficient[1] - (1.96* as.numeric(summary(exp_6a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[1] <- plot.data$Coefficient[1] + (1.96* as.numeric(summary(exp_6a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[2] <- as.numeric(summary(exp_7a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[2] <- plot.data$Coefficient[2] - (1.96* as.numeric(summary(exp_7a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[2] <- plot.data$Coefficient[2] + (1.96* as.numeric(summary(exp_7a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[3] <- as.numeric(summary(exp_8a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[3] <- plot.data$Coefficient[3] - (1.96* as.numeric(summary(exp_8a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[3] <- plot.data$Coefficient[3] + (1.96* as.numeric(summary(exp_8a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[4] <- as.numeric(summary(exp_1a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[4] <- plot.data$Coefficient[4] - (1.96* as.numeric(summary(exp_1a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[4] <- plot.data$Coefficient[4] + (1.96* as.numeric(summary(exp_1a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[5] <- as.numeric(summary(exp_4a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[5] <- plot.data$Coefficient[5] - (1.96* as.numeric(summary(exp_4a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[5] <- plot.data$Coefficient[5] + (1.96* as.numeric(summary(exp_4a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[6] <- as.numeric(summary(exp_5a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[6] <- plot.data$Coefficient[6] - (1.96* as.numeric(summary(exp_5a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[6] <- plot.data$Coefficient[6] + (1.96* as.numeric(summary(exp_5a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[7] <- as.numeric(summary(exp_2a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[7] <- plot.data$Coefficient[7] - (1.96* as.numeric(summary(exp_2a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[7] <- plot.data$Coefficient[7] + (1.96* as.numeric(summary(exp_2a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Coefficient[8] <- as.numeric(summary(exp_3a)$coefficients[c("postwar"),c("Estimate")])
plot.data$Lower[8] <- plot.data$Coefficient[8] - (1.96* as.numeric(summary(exp_3a)$coefficients[c("postwar"),c("Std. Error")]))
plot.data$Upper[8] <- plot.data$Coefficient[8] + (1.96* as.numeric(summary(exp_3a)$coefficients[c("postwar"),c("Std. Error")]))

plot.data$Lower <- as.numeric(plot.data$Lower)
plot.data$Upper <- as.numeric(plot.data$Upper)

plot.data$Variable <- as.factor(plot.data$Variable)
plot.data$Variable <- reorder(plot.data$Variable, plot.data$Coefficient)

table1_figure <- 
ggplot(data=plot.data, aes(x=Variable, y=Coefficient)) +
  geom_hline(yintercept=0, linetype="dotted") +
  geom_segment(data=plot.data, aes(x=Variable, xend=Variable, y=Lower, yend=Upper)) +
  geom_point(size=3) +
  geom_point(size=1, color="white") +
  coord_flip(ylim=c(-0.6,0.6)) +
  theme_classic() +
  scale_y_continuous(labels=function(x) paste0(x*100,"%" ) ) +
  xlab("") +
  ylab("Pre/Post Invasion Estimated Shift")

# save.plot("Table1_Figure.png", width=4.5, height=3)

save.tiff(table1_figure, "Table1_Figure.tiff",
          w=4.5, h=3)


## Pre and Post Models -- Robustness Check (1)
## exclude up to feb 20th (end of winter Olympics)

exp_13a <-  lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")), china_milpower ~ postwar + factor(country), weights=weight )
exp_13b <-  lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")), china_milpower ~ postwar + age + age_square + female + factor(country) + education + income  , weights=weight )

exp_14a <-  lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")), fav_china ~ postwar + factor(country) , weights=weight )
exp_14b <-  lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")), fav_china ~ postwar + age + age_square + female + factor(country) + education + income , weights=weight )
exp_14c <-  lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")), fav_china ~ postwar + age + age_square + female + factor(country) + education + income + fav_baseline, weights=weight )

exp_15c <- lm(data=subset(my.data, sample==1  & date>as.Date("2022-02-20")),  confid_xi ~ postwar + age + age_square + female + factor(country) + education + income  + confid_baseline 
              , weights=weight )
exp_15b <- lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")),  confid_xi ~ postwar + age + age_square + female + factor(country) + education + income 
              , weights=weight )
exp_15a <- lm(data=subset(my.data, sample==1 & date>as.Date("2022-02-20")),  confid_xi ~ postwar + factor(country) , weights=weight )

# get this to LaTeX
stargazer(exp_13a, exp_13b, exp_14a, exp_14b, exp_14c)  




## sample descriptives - pre and post 

table(data2022$sample)

data2022$fav_us_01 <- NA 
data2022$fav_us_01[data2022$fav_us==1] <- 0 
data2022$fav_us_01[data2022$fav_us==2] <- 0 
data2022$fav_us_01[data2022$fav_us==3] <- 1 
data2022$fav_us_01[data2022$fav_us==4] <- 1 


data2022$fav_ru_01 <- NA 
data2022$fav_ru_01[data2022$fav_russia==1] <- 0 
data2022$fav_ru_01[data2022$fav_russia==2] <- 0 
data2022$fav_ru_01[data2022$fav_russia==3] <- 1 
data2022$fav_ru_01[data2022$fav_russia==4] <- 1

data2022$confid_putin_01 <- NA 
data2022$confid_putin_01[data2022$confid_putin==1] <- 0 
data2022$confid_putin_01[data2022$confid_putin==2] <- 0 
data2022$confid_putin_01[data2022$confid_putin==3] <- 1 
data2022$confid_putin_01[data2022$confid_putin==4] <- 1

data2022$elderly <- NA 
data2022$elderly[data2022$age>70] <- 1 
data2022$elderly[data2022$age<71] <- 0 

data2022$income_high <- NA 
data2022$income_high[data2022$income>3] <- 1 
data2022$income_high[data2022$income<4] <- 0 

data2022$university <- NA 
data2022$university[data2022$education==3] <- 1 
data2022$university[data2022$education<3] <- 0 

data2022$demsat01 <- NA 
data2022$demsat01[data2022$demsat<3] <- 0 
data2022$demsat01[data2022$demsat>2] <- 1 

data2022$leftwing <- NA 
data2022$leftwing[data2022$ideology<4] <- 1 
data2022$leftwing[data2022$ideology>3] <- 0 


# check codings correct for recoded/added variables.

with(data2022, cor(demsat, demsat01, use="complete.obs") ) # ok
with(data2022, cor(fav_russia, confid_putin_01, use="complete.obs") ) # ok
with(data2022, cor(confid_putin, confid_putin_01, use="complete.obs") ) # ok
with(data2022, cor(fav_ru_01, fav_russia, use="complete.obs") ) # ok
with(data2022, cor(fav_ru_01, confid_putin, use="complete.obs") ) # ok
with(data2022, cor(fav_us_01, fav_us, use="complete.obs") ) # ok
with(data2022, cor(fav_us_01, confid_biden, use="complete.obs") ) # ok


return_stats <- function(variable)
{
mu.post <-  mean(
    as.vector(unlist(subset(data2022, postwar==1 & sample==1)[,c(variable)]))
    , na.rm=T)
mu.pre <-  mean(
  as.vector(unlist(subset(data2022, postwar==0 & sample==1)[,c(variable)]))
  , na.rm=T)
sd.post <- sd(
  as.vector(unlist(subset(data2022, postwar==1 & sample==1)[,c(variable)]))
  , na.rm=T)
sd.pre <- sd(
  as.vector(unlist(subset(data2022, postwar==0 & sample==1)[,c(variable)]))
  , na.rm=T)
count.post <- count(
  !is.na(unlist(subset(data2022, postwar==1 & sample==1)[,c(variable )]) ))[2,2]
count.pre <- count(
  !is.na(unlist(subset(data2022, postwar==0 & sample==1)[,c(variable )]) ))[2,2]
se.pre <- sd.pre / sqrt(count.pre)
se.post <- sd.post / sqrt(count.post)

foo <- as.data.frame(cbind(mu.pre, mu.pre - (se.pre*1.96), mu.pre + (se.pre*1.96), 
             mu.post, mu.post-se.post*1.96, mu.post+se.post*1.96))

names(foo)[1] <- paste("pre","mu",sep="_")
names(foo)[2] <- paste("pre", "lower",sep="_")
names(foo)[3] <- paste("pre", "upper",sep="_")
names(foo)[4] <- paste("post", "mu",sep="_")
names(foo)[5] <- paste("post", "lower",sep="_")
names(foo)[6] <- paste("post", "upper",sep="_")

rownames(foo) <- c(variable)

return( foo )
}



foo <- 

rbind(
return_stats("income_high"),
return_stats("female"),
return_stats("university"),
return_stats("elderly"),
return_stats("leftwing"),
return_stats("fav_ru_01"),
return_stats("confid_putin_01")
)

foo$vars <- rownames(foo)

foo$dv <- 0 
foo$dv[rownames(foo)=="fav_ru_01"] <- 1
foo$dv[rownames(foo)=="confid_putin_01"] <- 1


foo$vars <- c("Income (Upper Quintiles)", "Sex (Female)", "Education (University)",
              "Age (Elderly)","Ideology (Left/Center)",
              "Favorable View of Russia", "Confidence in Vladimir Putin")

foo$vars <- as.factor(foo$vars)
foo$vars <- reorder(foo$vars, -foo$dv)

descriptives_prepost <- 
ggplot(data=foo, aes(x=vars, y=pre_mu) ) + 
   geom_point(color="grey", size=3) +
  geom_segment(data=foo, aes(x=vars,xend=vars, y=post_lower, yend=post_upper),
               color="black", size=1 ) + 
  geom_point(data=foo, aes(x=vars, y=post_mu), color="black", size=3 ) + 
  geom_point(data=foo, aes(x=vars, y=post_mu), color="white", size=1.5 ) + 
  geom_segment(data=foo, aes(x=vars,xend=vars, y=pre_lower, yend=pre_upper),color="grey" ) + 
  coord_flip() +
   theme_classic() +
   scale_y_continuous(limits=c(0,1.1), breaks=c(0,0.25,0.5,0.75,1), 
                      labels=c("0%","25%","50%","75%","100%")) +
  xlab("") + ylab("") + 
  annotate(geom="rect", ymin=0.65, ymax=1.1, xmin=0.9, xmax=2.1, fill="grey", color="#EEEEEE", alpha=0.15, linewidth=0.35) +
  annotate(geom="point", y=0.7, x=1.25, size=3) + 
  annotate(geom="point", y=0.7, x=1.25, size=1.5, color="white") + 
  annotate(geom="point", y=0.7, x=1.75, size=3, color="grey") + 
  annotate(geom="text", y=0.75, x=1.25, size=3, label="Post-Invasion", hjust=0) +
  annotate(geom="text", y=0.75, x=1.75, size=3, label="Pre-Invasion", hjust=0)

# save.plot("descriptives_prepost.png", width=4.5, height=5)

save.tiff(descriptives_prepost, "descriptives_prepost.tiff",
          w=4.5, h=5)


## check effect if Sweden removed from sample? 

data2022$sample[data2022$country=="Sweden"] <- 0


foo <- 
  
  rbind(
    return_stats("income_high"),
    return_stats("female"),
    return_stats("university"),
    return_stats("elderly"),
    return_stats("leftwing"),
    return_stats("fav_ru_01"),
    return_stats("confid_putin_01")
  )

foo$vars <- rownames(foo)

foo$dv <- 0 
foo$dv[rownames(foo)=="fav_ru_01"] <- 1
foo$dv[rownames(foo)=="confid_putin_01"] <- 1


foo$vars <- c("Income (Upper Quintiles)", "Sex (Female)", "Education (University)",
              "Age (Elderly)","Ideology (Left/Center)",
              "Favorable View of Russia", "Confidence in Vladimir Putin")

foo$vars <- as.factor(foo$vars)
foo$vars <- reorder(foo$vars, -foo$dv)

ggplot(data=foo, aes(x=vars, y=pre_mu) ) + 
  geom_point(color="grey", size=3) +
  geom_segment(data=foo, aes(x=vars,xend=vars, y=post_lower, yend=post_upper),
               color="black", size=1 ) + 
  geom_point(data=foo, aes(x=vars, y=post_mu), color="black", size=3 ) + 
  geom_point(data=foo, aes(x=vars, y=post_mu), color="white", size=1.5 ) + 
  geom_segment(data=foo, aes(x=vars,xend=vars, y=pre_lower, yend=pre_upper),color="grey" ) + 
  coord_flip() +
  theme_classic() +
  scale_y_continuous(limits=c(0,1), breaks=c(0,0.25,0.5,0.75,1), 
                     labels=c("0%","25%","50%","75%","100%")) +
  xlab("") + ylab("")

## this figure did not make final cut:
## but good to include code in case anyone wanted to check. 


## not much difference, but best to check in order to avoid complaints.


ttest_table <-
as.data.frame(table(subset(data2022, sample==1)$country) )

names(ttest_table)[1] <- "country"

ttest_table$t_score <- NA 
ttest_table$p_value <- NA 


for(i in 1:length(ttest_table$country))
{
  ttest_table$t_score[i] <- 
  t.test(fav_ru_01 ~ postwar, data = subset(data2022, country==ttest_table$country[i]) )$stat

  ttest_table$p_value[i] <- 
    t.test(fav_ru_01 ~ postwar, data = subset(data2022, country==ttest_table$country[i]) )$p.value
  
}

ttest_table$lower <- NA 
ttest_table$upper <- NA 

for(i in 1:length(ttest_table$country))
{
  ttest_table$lower[i] <- 
    t.test(fav_ru_01 ~ postwar, data = subset(data2022, country==ttest_table$country[i]) )$conf.int[1]
  ttest_table$upper[i] <- 
    t.test(fav_ru_01 ~ postwar, data = subset(data2022, country==ttest_table$country[i]) )$conf.int[2]
}
  
ttest_table$p_value <- round(ttest_table$p_value, 4)

# table we reported in the appendix 
ttest_table


## table of descriptive statistics for all variables used 


return_descriptives <- function(variable)
{
  mu <-  mean(as.vector(unlist(pewdata[,c(variable)])), na.rm=T)
  sd <- sd(as.vector(unlist(pewdata[,c(variable)])), na.rm=T)
  count <- count(!is.na(unlist(pewdata[,c(variable )]) ))[2,2]
  se <- sd / sqrt(count-1)  
  ci.low <- mu - 1.96*se
  ci.upp <- mu + 1.96*se
  foo <- as.data.frame(cbind(mu, sd, count, ci.low, ci.upp))
  rownames(foo)[1] <- variable
  
  return(foo)
}  


## I probably want to convert to binary variables here for pewdata too.  

pewdata$fav_us01 <- ifelse(pewdata$fav_us>2, 1, 0)
pewdata$fav_china01 <- ifelse(pewdata$fav_china>2, 1, 0)
pewdata$fav_russia01 <- ifelse(pewdata$fav_russia>2, 1, 0)
pewdata$fav_eu01 <- ifelse(pewdata$fav_eu>2, 1, 0)
pewdata$university <- ifelse(pewdata$education==3, 1, 0)

pewdata$demsat01 <- ifelse(pewdata$demsat>2, 1, 0)
pewdata$religion_import01 <- ifelse(pewdata$religion_import>2, 1, 0)
pewdata$econ_sit01 <- ifelse(pewdata$econ_sit<3, 1, 0)

# check polarities and recodings  
# we already confirmed that fav_* were correctly scaled so basic check here
with(pewdata, cor(demsat, demsat01, use="complete.obs") ) # ok
with(pewdata, cor(fav_us01, fav_us, use="complete.obs") ) # ok
with(pewdata, cor(fav_china01, fav_china, use="complete.obs") ) # ok
with(pewdata, cor(fav_russia01, fav_russia, use="complete.obs") ) # ok
with(pewdata, cor(fav_eu01, fav_eu, use="complete.obs") ) # ok

# check scaled low to high 
with(pewdata, table(country, religion_import) ) # ok
with(pewdata, cor(religion_import, religion_import01, use="complete.obs") ) # ok
# check scaling. For econ_sit it should be high to low 
with(pewdata, table(country, econ_sit) ) # ok
# then we inverted this when generating econ_sit01 (0 = bad, 1= good)
with(pewdata, cor(econ_sit01, econ_sit, use="complete.obs") ) # so negative = ok


foo <- rbind(
  return_descriptives("demsat01"),
  return_descriptives("fav_us01"),
  return_descriptives("fav_china01"),
  return_descriptives("fav_russia01"),
  return_descriptives("fav_eu01"),
  return_descriptives("age"),
  return_descriptives("female"),
  return_descriptives("university"),
  return_descriptives("ideology"),
  return_descriptives("income"),
  return_descriptives("religion_import01"),
  return_descriptives("econ_sit01")
)


# descriptives table 
foo




## Load 2023 data separately.  
## We needed a new version of the file to have the support for democracy variable. 
## This is the 2023 file only, with all 2023 variables. 


pewdata2023 <- read_dta("pew2023data_precode.dta")

# once again have to check the scales/codings are all correct for this version of the file:  
# also checked via the variable labels, but best to be sure in case of any precoding in Stata that didn't update the labels.

# Pew by default code a lot of variables positive to negative, so almost all vars need rescaling. 
# when we load from scratch the raw download for the latest (2023) dataset.

with(pewdata2023, table(country, demsat) ) # needs rescaling 
pewdata2023$demsat <- 5 - pewdata2023$demsat 

with(pewdata2023, table(country, polsys_republic) ) # needs rescaling 
pewdata2023$polsys_republic <- 5 - pewdata2023$polsys_republic 

with(pewdata2023, cor(demsat, polsys_republic, use="complete.obs") ) # ok

with(pewdata2023, table(country, fav_us) ) # needs rescaling 
pewdata2023$fav_us <- 5 - pewdata2023$fav_us 
pewdata2023$fav_china <- 5 - pewdata2023$fav_china 
pewdata2023$fav_eu <- 5 - pewdata2023$fav_eu
with(pewdata2023, cor(demsat, fav_eu, use="complete.obs") ) # ok
with(pewdata2023, cor(fav_us, fav_eu, use="complete.obs") ) # ok
pewdata2023$fav_un <- 5 - pewdata2023$fav_un
with(pewdata2023, cor(fav_un, fav_eu, use="complete.obs") ) # ok
pewdata2023$fav_nato <- 5 - pewdata2023$fav_nato
with(pewdata2023, cor(fav_us, fav_nato, use="complete.obs") ) # ok
pewdata2023$fav_russia <- 5 - pewdata2023$fav_russia
pewdata2023$confid_biden <- 5 - pewdata2023$confid_biden 
with(pewdata2023, cor(fav_us, confid_biden, use="complete.obs") ) # ok
pewdata2023$confid_xi <- 5 - pewdata2023$confid_xi
with(pewdata2023, cor(fav_china, confid_xi, use="complete.obs") ) # ok
pewdata2023$confid_putin <- 5 - pewdata2023$confid_putin
with(pewdata2023, cor(fav_russia, confid_putin, use="complete.obs") ) # ok
pewdata2023$confid_macron <- 5 - pewdata2023$confid_macron
with(pewdata2023, cor(fav_eu, confid_macron, use="complete.obs") ) # ok
pewdata2023$confid_scholz <- 5 - pewdata2023$confid_scholz
with(pewdata2023, cor(fav_eu, confid_scholz, use="complete.obs") ) # ok
pewdata2023$confid_zelenskyy <- 5 - pewdata2023$confid_zelenskyy
with(pewdata2023, cor(fav_nato, confid_zelenskyy, use="complete.obs") ) # ok
with(pewdata2023, cor(confid_putin, confid_zelenskyy, use="complete.obs") ) # negative as expected
pewdata2023$polsys_directdemocracy <- 5 - pewdata2023$polsys_directdemocracy 
with(pewdata2023, cor(polsys_directdemocracy, polsys_republic, use="complete.obs") ) # ok
# no reason to think this is miscoded but can check. Sex coded as female = 2
table(pewdata2023$sex,pewdata2023$female) # ok

with(pewdata2023, table(country, abortion_legal) ) # needs rescaling 
pewdata2023$abortion_legal <- 5 - pewdata2023$abortion_legal 

# for visualisation, better to have a dummy age variable, as otherwise coefficient misleadingly small 
# as a report of age effects

table(pewdata2023$age) # checking 98 and 99 values already removed as NAs 

pewdata2023$older <- NA 
pewdata2023$older[pewdata2023$age>35 & !is.na(pewdata2023$age)] <- 1
pewdata2023$older[pewdata2023$age<36 & !is.na(pewdata2023$age)] <- 0

table(pewdata2023$age,pewdata2023$older) # ok 


return_std_coef <- function(variable="demsat")
{
foo <- summary(lm(data=pewdata2023, scale(polsys_republic) ~ scale(  eval(parse(text = paste(variable))) ) + factor(country) ) )
coef <- foo$coefficients[2,1]
se <- foo$coefficients[2,2]
bar <- as.data.frame(cbind(coef,se) )
rownames(bar) <- variable
return(bar) 
}

democracy_coefs <- 
  rbind(return_std_coef("demsat"),
return_std_coef("fav_us"),
return_std_coef("fav_china"),
return_std_coef("fav_eu"), 
return_std_coef("fav_un"),
return_std_coef("fav_nato"),
return_std_coef("fav_russia"),
return_std_coef("confid_biden"),
return_std_coef("confid_xi"),
return_std_coef("confid_putin"),
return_std_coef("confid_macron"),
#return_std_coef("confid_scholz"), # omitted for parsimony but leaving the code if anyone wants to check
return_std_coef("confid_zelenskyy"),
return_std_coef("polsys_directdemocracy"),
return_std_coef("female"),
return_std_coef("ideology"),
return_std_coef("abortion_legal"),
return_std_coef("older") )

democracy_coefs$lower <- democracy_coefs$coef - 1.96*democracy_coefs$se
democracy_coefs$upper <- democracy_coefs$coef + 1.96*democracy_coefs$se

democracy_coefs$name <- c("Satisfied with Democracy", "Favor the United States", "Favor China", 
                          "Favor the European Union", "Favor the United Nations", "Favor the NATO Alliance", "Favor Russia",
                          "Confidence in President Biden", "Confidence in President Xi", "Confidence in President Putin",
                          "Confidence in President Macron", "Confidence in President Zelenskyy",
                          "Approve of Direct Democracy", "Gender (Female = 1)", "Ideology (Left-Right)",
                          "Support Abortion Rights", "Age (Over 35 = 1)")


democracy_coefs$name <- factor(democracy_coefs$name)
democracy_coefs$name <- reorder(democracy_coefs$name, democracy_coefs$coef)

democracy_coefs_check <- 
ggplot(data=democracy_coefs, aes(x=name, y=coef)) + 
  geom_segment(data=democracy_coefs, aes(x=name, xend=name, y=lower, yend=upper)) + 
   geom_point(size=2.5) + 
  geom_point(size=1, color="white") + 
  coord_flip() + 
   xlab("") + ylab("Standardized Coefficient, Controlling for CFEs") +  
  theme_classic()

# save.plot("democracy_coefs_check.png", width=6, height=4.5)

save.tiff(democracy_coefs_check, "democracy_coefs_check.tiff",
           w=6, h=4.5)


## liberal democracies only robustness check. 

return_std_coef_d <- function(variable="demsat")
{
  foo <- summary(lm(data=subset(pewdata2023, country!="Hungary" & country!="India" & country!="Indonesia"
                                & country !="Kenya" & country!="Mexico" & country!="Nigeria" ),
          scale(polsys_republic) ~ scale(  eval(parse(text = paste(variable))) ) + factor(country) ) )
  coef <- foo$coefficients[2,1]
  se <- foo$coefficients[2,2]
  bar <- as.data.frame(cbind(coef,se) )
  rownames(bar) <- variable
  return(bar) 
}

democracy_coefs <- 
  rbind(return_std_coef_d("demsat"),
        return_std_coef_d("fav_us"),
        return_std_coef_d("fav_china"),
        return_std_coef_d("fav_eu"), 
        return_std_coef_d("fav_un"),
        return_std_coef_d("fav_nato"),
        return_std_coef_d("fav_russia"),
        return_std_coef_d("confid_biden"),
        return_std_coef_d("confid_xi"),
        return_std_coef_d("confid_putin"),
        return_std_coef_d("confid_macron"),
        return_std_coef_d("confid_scholz"),  # included here as just a robustness check 
        return_std_coef_d("confid_zelenskyy"),
        return_std_coef_d("polsys_directdemocracy"),
        return_std_coef_d("female"),
        return_std_coef_d("ideology"),
        return_std_coef_d("abortion_legal"),
        return_std_coef_d("older") )

democracy_coefs$lower <- democracy_coefs$coef - 1.96*democracy_coefs$se
democracy_coefs$upper <- democracy_coefs$coef + 1.96*democracy_coefs$se

democracy_coefs$name <- c("Satisfied with Democracy", "Favor the United States", "Favor China", 
                          "Favor the European Union", "Favor the United Nations", "Favor the NATO Alliance", "Favor Russia",
                          "Confidence in President Biden", "Confidence in President Xi", "Confidence in President Putin",
                          "Confidence in President Macron", "Confidence in Chancellor Scholz", "Confidence in President Zelenskyy",
                          "Approve of Direct Democracy", "Gender (Female = 1)", "Ideology (Left-Right)",
                          "Support Abortion Rights", "Age (Over 35 = 1)")


democracy_coefs$name <- factor(democracy_coefs$name)
democracy_coefs$name <- reorder(democracy_coefs$name, democracy_coefs$coef)

democracy_coefs_check_demonly <- 
ggplot(data=democracy_coefs, aes(x=name, y=coef)) + 
  geom_segment(data=democracy_coefs, aes(x=name, xend=name, y=lower, yend=upper)) + 
  geom_point(size=2.5) + 
  geom_point(size=1, color="white") + 
  coord_flip() + 
  xlab("") + ylab("Standardized Coefficient, Controlling for CFEs") +  
  theme_classic()

# save.plot("democracy_coefs_check_demonly.png", width=6, height=4.5)

save.tiff(democracy_coefs_check_demonly, "democracy_coefs_check_demonly.tiff",
          w=6, h=4.5)



## support for democracy coefs chart. 
## We exported the coefficients already to a csv file to make this simpler. 

demsup_coefs <- read.csv("demsup_2017_2023_coefs.csv")

demsup_coefs$lower <- demsup_coefs$coef - (1.645*demsup_coefs$se)
demsup_coefs$upper <- demsup_coefs$coef + (1.645*demsup_coefs$se)

demsup_coefs$country <- c("United States", "China", "Russia", "European Union",
                          "United States", "China", "Russia", "European Union" )
  
demsup_coefs$country <- factor(demsup_coefs$country )
demsup_coefs$country <- reorder(demsup_coefs$country , demsup_coefs$coef)
  
demsup_coefs$end <- NA 
demsup_coefs$end[demsup_coefs$year==2017] <- demsup_coefs$coef[demsup_coefs$year==2023]

demsup_2017_2023_check <- 
ggplot(data=demsup_coefs, aes(x=country, y=coef)) + 
  geom_segment(data=subset(demsup_coefs, year==2023), 
               aes(x=country, xend=country, y=lower, yend=upper), color="black", size=0.5) + 
  geom_segment(data=subset(demsup_coefs, year == 2017 & end>coef),
               aes(x=country, xend=country, y=coef, yend=coef + 0.025), size = 1.75, alpha=1,
               color="darkgrey", arrow = arrow(type = "closed",length = unit(10, "points")) ) +
  geom_segment(data=subset(demsup_coefs, year == 2017 & end<coef),
               aes(x=country, xend=country, y=coef, yend=end + 0.005), size = 1.75, alpha=1,
               color="darkgrey", arrow = arrow(type = "closed",length = unit(10, "points")) ) +
  geom_point(size=5.5) + 
  geom_point(data=subset(demsup_coefs, year==2017),
             aes(x=country, y=coef), color="darkgrey", size=5.5) + 
  geom_point(size=1, color="white") + 
  coord_flip() +
  xlab("") + ylab("Standardized Coefficient, Controlling for CFEs") +  
  theme_classic() + 
  theme(text = element_text(size = 15) ) 


# save.plot("demsup_2017_2023_check.png", width=10, height=3)

save.tiff(demsup_2017_2023_check, "demsup_2017_2023_check.tiff",
          w=10, h=3)




## MLM Models Chart for Country Leader Favorability Over Time


ml_14_1 <- lmer(data=subset(pewdata, constant==1), confid_putin ~  age + age_square + female + factor(country) + factor(year) + education + income + religion_import +  econ_sit + confid_baseline +
                  ( 1 + demsat | year ), weights=weight )


poo1 <- ranef(ml_14_1)$year
poo1 <- adf(poo1)
poo1$year <- as.numeric( row.names(poo1) )

poo1$se.demsat <- se.ranef(ml_14_1)$year[,"demsat"]
poo1$demsat.lower <- poo1$demsat - 1.645*poo1$se.demsat
poo1$demsat.upper <- poo1$demsat + 1.645*poo1$se.demsat


ml_14_2 <- lmer(data=subset(pewdata, constant==1), confid_xi ~  age + age_square + female + factor(country) + factor(year) + education + income + religion_import +  econ_sit + confid_baseline +
                  ( 1 + demsat | year ), weights=weight )


xoo1 <- ranef(ml_14_2)$year
xoo1 <- adf(xoo1)
xoo1$year <- as.numeric( row.names(xoo1) )

xoo1$se.demsat <- se.ranef(ml_14_2)$year[,"demsat"]
xoo1$demsat.lower <- xoo1$demsat - 1.645*xoo1$se.demsat
xoo1$demsat.upper <- xoo1$demsat + 1.645*xoo1$se.demsat


pewdata$confid_uspres <- rowMeans( cbind(as.numeric(pewdata$confid_trump), as.numeric(pewdata$confid_biden), as.numeric(pewdata$confid_obama)), na.rm=T )


ml_14_3 <- lmer(data=subset(pewdata, constant==1), confid_uspres ~  age + age_square + female + factor(country) + factor(year) + education + income + religion_import +  econ_sit + confid_baseline +
                  ( 1 + demsat | year ), weights=weight )


too1 <- ranef(ml_14_3)$year
too1 <- adf(too1)
too1$year <- as.numeric( row.names(too1) )

too1$se.demsat <- se.ranef(ml_14_3)$year[,"demsat"]
too1$demsat.lower <- too1$demsat - 1.645*too1$se.demsat
too1$demsat.upper <- too1$demsat + 1.645*too1$se.demsat

MLM_slopes_demsat_bycountry_leaders_new_ccsample <-
ggplot(subset(poo1, year>2017 ), aes(x=year, y=demsat)) +
  geom_hline(yintercept=0, size=0.5, color="darkgrey") + 
  geom_ribbon(aes(ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_ribbon(data=subset(xoo1,year>2017 ), aes(x=year, ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_ribbon(data=subset(too1,year>2017 ), aes(x=year, ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) +
  geom_segment(data= subset(poo1, year==2023 | year==2018), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data= subset(xoo1, year==2023 | year==2018), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data= subset(too1, year==2023 | year==2018), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_line(color="darkgrey", size=2) + 
  geom_line(data=subset(xoo1, year>2017 ), aes(year, demsat), size=2) + 
  geom_line(data=subset(too1, year>2017 ), aes(year, demsat  ), size=2, color="#444444") +
  geom_point(data=subset(poo1, year==2023 | year==2018), aes(year,demsat)) + 
  geom_point(data=subset(xoo1, year==2023 | year==2018), aes(year,demsat)) + 
  geom_point(data=subset(too1, year==2023 | year==2018), aes(year,demsat)) + 
  xlab("") + 
  ylab("Satisfaction with Democracy Estimated Effect") + 
  theme_article(base_size=15) +  
  theme(panel.grid.minor.y = element_blank()) + 
  scale_x_continuous( breaks=c(2010, 2012,2014,2016,2018,2020,2022), limits=c(2012,2024)) + 
  scale_y_continuous( breaks=c(-0.2,-0.1,0,0.1,0.2,0.3,0.4), labels=c("-0.20","-0.10","0.00","+0.10","+0.20","+0.30","+0.40")) + 
  coord_cartesian(xlim=c(2017.5,2024.5)) + 
  annotate(geom="text", label="President Xi", hjust=0, x=2023.2, y= tail(xoo1$demsat,1) ) +
  annotate(geom="text", label="President Putin", hjust=0, x=2023.2, y= tail(poo1$demsat,1)  ) +
  annotate(geom="text", label="US President", hjust=0, x=2023.2, y= tail(too1$demsat,1)  ) 


# save.plot("MLM_slopes_demsat_bycountry_leaders_new_ccsample.png", width=9)

save.tiff(MLM_slopes_demsat_bycountry_leaders_new_ccsample, 
          "MLM_slopes_demsat_bycountry_leaders_new_ccsample.tiff",
          w=9, h=6)



## let's check religion over time by country opinion.


ml_9_5 <- lmer(data=pewdata , prowestern ~ demsat + age + age_square + female + factor(country)  + education + income  +  econ_sit  + fav_baseline +
                 ( 1 + religion_import | year ), weights=weight )
 
roo5 <- ranef(ml_9_5)$year
roo5 <- adf(roo5)
roo5$year <- as.numeric( row.names(roo5) )

roo5$se.religion_import <- se.ranef(ml_9_5)$year[,"religion_import"]
roo5$religion_import.lower <- roo5$religion_import - 1.645*roo5$se.religion_import
roo5$religion_import.upper <- roo5$religion_import + 1.645*roo5$se.religion_import

Relig_prowestern <-
ggplot(subset(roo5, year>2016  ), aes(x=year, y=religion_import)) +
  geom_ribbon(aes(ymax=religion_import.upper, ymin=religion_import.lower ),alpha=.05 ) + 
  geom_line() + 
  geom_segment(data= subset(roo5, year== 2017 | year==2023), aes(x=year, xend=year, y=religion_import.lower, yend=religion_import.upper), color="darkgray" ) +
  geom_point(data=subset(roo5, year==2023 | year==2017), aes(year,religion_import)) + 
  theme_article(base_size=15) +  
  xlab("") + 
  coord_cartesian(ylim=c(-0.2, 0.15), xlim=c(2016.5,2025)) + 
  ylab("Religiosity Estimated Effect") +
  scale_x_continuous( breaks=c(2010, 2012,2014,2016,2018,2020,2022), limits=c(2012,2024)) + 
  scale_y_continuous( breaks=c(-0.2,-0.1,0,0.1,0.2,0.3,0.4), labels=c("-0.20","-0.10","0.00","+0.10","+0.20","+0.30","+0.40")) + 
  annotate(geom="text", label="Pro-Western", hjust=0, x=2023.2, y= tail(roo5$religion_import,1) )


# save.plot("Relig_prowestern.png", width=9)

save.tiff(Relig_prowestern, 
          "Relig_prowestern.tiff",
          w=9, h=6)


### Estimating Separately for the Global South 

ml_10_1 <- lmer(data= subset(my.data, globalsouth==1), prowestern ~   age + age_square + female  + factor(country) + education + income + religion_import + econ_sit + fav_baseline +
                  ( 1 + demsat | year ), weights=weight )


zoo1 <- ranef(ml_10_1)$year
zoo1 <- adf(zoo1)
zoo1$year <- as.numeric( row.names(zoo1) )

zoo1$se.demsat <- se.ranef(ml_10_1)$year[,"demsat"]
zoo1$demsat.lower <- zoo1$demsat - 1.645*zoo1$se.demsat
zoo1$demsat.upper <- zoo1$demsat + 1.645*zoo1$se.demsat


ml_10_2 <- lmer(data= subset(my.data, globalsouth==0), prowestern ~   age + age_square + female  + factor(country) + education + income + religion_import + econ_sit + fav_baseline +
                  ( 1 + demsat | year ), weights=weight )


zoo2 <- ranef(ml_10_2)$year
zoo2 <- adf(zoo2)
zoo2$year <- as.numeric( row.names(zoo2) )

zoo2$se.demsat <- se.ranef(ml_10_2)$year[,"demsat"]
zoo2$demsat.lower <- zoo2$demsat - 1.645*zoo2$se.demsat
zoo2$demsat.upper <- zoo2$demsat + 1.645*zoo2$se.demsat
zoo2 <- subset(zoo2, year!=2018)

Global_North_South_prowestern <- 
ggplot(subset(zoo1, year>2012 ), aes(x=year, y=demsat)) +
  geom_ribbon(aes(ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) + 
  geom_line() + 
  geom_ribbon(data=subset(zoo2,year>2016),aes(ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) + 
  geom_line(data=subset(zoo2,year>2016 ), aes(x=year,y=demsat)) + 
  geom_segment(data= subset(zoo1, year== 2013 | year==2022), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data= subset(zoo2, year== 2017 | year==2022), aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_point(data=subset(zoo1, year>2012), aes(year,demsat)) + 
  geom_point(data=subset(zoo2, year>2016), aes(year,demsat)) + 
  theme_article(base_size=15) +  
  xlab("") + 
  coord_cartesian(ylim=c(-0.2, 0.5), xlim=c(2012.5,2024)) + 
  ylab("Satisfaction with Democracy Estimated Effect") +
  scale_x_continuous( breaks=c(2010, 2012,2014,2016,2018,2020,2022), limits=c(2012,2023)) + 
  scale_y_continuous( breaks=c(-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5), labels=c("-0.20","-0.10","0.00","+0.10","+0.20","+0.30","+0.40","+0.50")) + 
  annotate(geom="text", label="Global\nSouth", hjust=0, x=2022.2, y= tail(zoo1$demsat,1) ) +
  annotate(geom="text", label="Global\nNorth", hjust=0, x=2022.2, y= tail(zoo2$demsat,1) )

# save.plot("Global_North_South_prowestern.png", width=9)

save.tiff(Global_North_South_prowestern, 
          "Global_North_South_prowestern.tiff",
          w=9, h=6)



##  favorability to NATO over time. 


ml_11_1 <- lmer(data= subset(my.data, ally=="NATO" &  constant==1), fav_nato ~   age + age_square + female  + factor(country) + education + income + religion_import + econ_sit + fav_baseline +
                  ( 1 + demsat | year ), weights=weight )

bung <- ranef(ml_11_1)$year
bung <- adf(bung)
bung$year <- as.numeric( row.names(bung) )

bung$se.demsat <- se.ranef(ml_11_1)$year[,"demsat"]
bung$demsat.lower <- bung$demsat - 1.645*bung$se.demsat
bung$demsat.upper <- bung$demsat + 1.645*bung$se.demsat


ml_11_2 <- lmer(data= subset(my.data, ally=="NATO" & constant==1), fav_un ~   age + age_square + female  + factor(country) + education + income + religion_import + econ_sit + fav_baseline +
                  ( 1 + demsat | year ), weights=weight )

sung <- ranef(ml_11_2)$year
sung <- adf(sung)
sung$year <- as.numeric( row.names(sung) )

sung$se.demsat <- se.ranef(ml_11_2)$year[,"demsat"]
sung$demsat.lower <- sung$demsat - 1.645*sung$se.demsat
sung$demsat.upper <- sung$demsat + 1.645*sung$se.demsat

mlm_nato <- 
ggplot(subset(bung, year>2016 ), aes(x=year, y=demsat)) +
  geom_ribbon(aes(ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) + 
  geom_ribbon(data=subset(sung, year>2016), aes(ymax=demsat.upper, ymin=demsat.lower ),alpha=.05 ) + 
  geom_line() + 
  geom_line(data=subset(sung, year>2016), aes(x=year, y=demsat)) + 
  geom_segment( aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_segment(data=subset(sung, year>2016),  aes(x=year, xend=year, y=demsat.lower, yend=demsat.upper), color="darkgray" ) +
  geom_point() + 
  geom_point(data=subset(sung, year>2016), aes(x=year, y=demsat)) + 
  theme_article(base_size=15) +  
  xlab("") + 
  coord_cartesian(ylim=c(0, 0.3), xlim=c(2016.5,2024.5)) + 
  ylab("Satisfaction with Democracy Estimated Effect") +
  scale_x_continuous( breaks=c(2010, 2012,2014,2016,2018,2020,2022), limits=c(2012,2024)) + 
  scale_y_continuous( breaks=c(-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5), labels=c("-0.20","-0.10","0.00","+0.10","+0.20","+0.30","+0.40","+0.50")) + 
  annotate(geom="text", label="Favorable to NATO", hjust=0, x=2022.2, y= tail(bung$demsat,1) ) + 
  annotate(geom="text", label="Favorable to the UN", hjust=0, x=2022.2, y= tail(sung$demsat,1) )

# save.plot("mlm_nato.png", width=9)

save.tiff(mlm_nato, 
          "mlm_nato.tiff",
          w=9, h=6)
